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The growth rates and chemical ordering of ferroelectric alloys are studied with kinetic Monte Carlo 
(KMC) simulations using an electrostatic model with long-range Coulomb interactions, as a function 
of temperature, chemical composition, and substrate orientation. Crystal growth is characterized 
by thermodynamic processes involving adsorption and evaporation, with solid-on-solid restrictions 
and excluding diffusion. A KMC algorithm is formulated to simulate this model efficiently in the 
presence of long-range interactions. Simulations were carried out on Ba(Mg 1/ / 3 Nb2/3)03 (BMN) type 
materials. Compared to the simple rocksalt ordered structures, ordered BMN grows only at very low 
temperatures and only under finely tuned conditions. For materials with tetravalent compositions, 
such as (1 — a;)Ba(Mg 1 y 3 Nb 2 /3)0 3 +a;BaZr03 (BMN-BZ), the model does not incorporate tetravalent 
ions at low-temperature, exhibiting a phase-separated ground state instead. At higher temperatures, 
tetravalent ions can be incorporated, but the resulting crystals show no chemical ordering in the 
absence of diffusive mechanisms. 



I. INTRODUCTION 

Ferroelectric crystals are known for their important 
technological applications such as high-pcrmitivity di- 
electrics, piezoelectric sensors, transducers, and me- 
chanical actuators i Recently, single-crystal relaxor per- 
ovskites such as Pb(Zn 1/3 Nb 2/ 3)0 3 -PbTi03 (PZN-PT) 
and Pb(Mg 1/3 Nb 2/3 )03-PbTi03 (PMN-PT) were syn- 
thesized and found to exhibit ultrahigh stain and very 
large piezoelectric constants^ The structure of alloys like 
PMN-PT can be viewed as a perovskitc AB0 3 framework 
(a cubic lattice for the ideal perovskite crystal), with Pb 
ions on the A-site and a solid solution of (Mg +2 , Nb +5 , 
Ti +4 ) ions on the B-sites, with average +4 B-site ionic 
charge. Of course this is an idealized picture, neglect- 
ing vacancies, impurities, local structural distortions, and 
partial chemical ordering on the B-sites. 

Partial B-sitc chemical ordering is a common feature 
of the high-piezoelectric solid solutions. While random 
B-site ordering is observed in isoclectronic solid solu- 
tions like Pb(Zri_ x Ti x )0 3 (PZT), non-isoelectronic B- 
site solid solutions A(BB B )0 3 , with B-site cations 
from group II, IV, and V, often exhibit compositionally- 
dependent B-site chemical ordering. At 1640°C, when 
the tetravalent composition x is increased in (1 — 
x)Ba(Mg 1/3 Nb 2/3 )0 3 + x BaZr0 3 (BMN-BZ), the fol- 
lowing sequence of B-site ordering is observed: [lll]i : 2 
order for x < 5%; then [lll]i : i order for 5% < x < 25%; 
and finally disorder for larger x& The [1 1 1] i : 2 nota- 
tion refers to x-ray observation of alternating (3(3(3' [111] 
stacking of B-sites, where (3 and (3' denote average scat- 
tering sites. For example, in BMN-BZ with x = 0, 
one can identify (3 with Nb and (3' with Mg. The 
[111] i : i notation refers to x-ray observation of rocksalt- 
like alternating (3(3' [111] stacking of B-cations. In this 
case, the assignment of the (3 and (3' sites has been de- 
bated, as discussed below in connection with the space- 
charge and random-site models^. Other Ba-based per- 
ovskites, e.g., (l-x) Ba(Mgi/ 3 Ta 2 / 3 )0 3 + x BaZr0 3 



(BMT-BZ)^ (l-x) Ba(Mg 1/3 Nb 2/3 )0 3 + x BaZr0 3 
(BMN-BZ) ji display a similar sequence of B-site or- 
der. On the other hand, for Pb-based systems, e.g., (1- 
x) Pb(Mgi /3 Ta 2/3 )0 3 + x PbZr0 3 (PMT-PZ), [111] 1:2 
order is not observed at x = 0; instead, annealing be- 
tween 1325°C and 1350°C results in [lll]i : i order all the 
way down to x = Other Pb-based perovskitcs , e.g., 
Pb(Mg 1/3 Nb 2/3 )0 3 (PMN),££ display similar B-site or- 
dering. 

Since their discovery, growing large single crystals has 
been a major research goal, but this effort has been 
largely unsupported by theory, because of the difficulty 
in modelling and simulating the non-equilibrium pro- 
cesses occurring in nucleation and crystal growth in such 
complex materials. In this paper, we use kinetic Monte 
Carlo^ simulations of a simple effective Hamiltonian to 
model the growth process of these ferroelectric crystals. 

Given the ionic character of these materials, it is not 
surprising that the inclusion of Coulomb interactions 
has been found to be crucial in describing their prop- 
erties. A simple, purely electrostatic model introduced 
by Bellaiche and Vandcrbilt (BV)i£ has had consider- 
able success in explaining the observed equilibrium B- 
site chemical ordering in many perovskite alloys. The 
BV model only considers Coulomb interactions between 
point charges (+2, +5, +4, etc., representing the differ- 
ent atomic species) that reside on the B-sites, which are 
constrained to lie on an ideal cubic sublattice. 

This electrostatic model is the starting point of our 
growth simulations. Simplified models based on Ising 
like effective Hamiltonians H c fi have been used to model 
growth in simpler systemsii*^. These models often have 
only short-range interactions. To adapt the electrostatic 
model of BV to study crystal growth, we consider a 
slab-geometry with periodic boundary conditions in two- 
dimensions only. The slab is viewed as being embedded 
in a liquid-phase melt, which is parametrized by a chem- 
ical potential difference with the solid bulk phase. In 
our simulations, a solid-on-solid (SOS) restriction is im- 



2 



posed, which requires that adsorption only occur onto 
empty lattice sites directly above an occupied site, so 
void formation is neglected. In keeping with the simplic- 
ity of the model, diffusion in the bulk and at the surface 
is also neglected. 

The non-equilibrium dynamics of the growth process is 
modeled using the kinetic Monte Carlo (KMC) method^. 
The KMC algorithm introduced by Bortz, Kalos, and 
Lcbowitz (BKL)S has been quite successful in simulat- 
ing crystal growth in Ising-like models with short-range 
interactions between adatoms. We generalize the algo- 
rithm to efficiently handle the long-range interactions in 
the BV model. 

In this paper we study the properties of this mini- 
mal paradigm of the growth process to determine if it 
can yield insights into the physics of the observed B-sitc 
chemical ordering. Section II presents our theoretical ap- 
proach. The model of BV is reviewed, and our adaptation 
of it for the growth modeling is discussed, including the 
special handling of electrostatic interactions during the 
growth process. Our generalization and modification of 
the KMC algorithm are then described which allows effi- 
cient treatment of the long-range Coulomb interactions. 
Section III presents the results of our growth simulations 
for A(BB')0 3 and A(BB'b")0 3 crystals. To help under- 
stand our growth results for the latter systems, where a 
(typically small) fraction of tetravalent B ions are mixed 
in, we also carry out total energy calculations to study 
their stability. Finally in section IV we further discuss 
our results and prospects for the model. In the appendix, 
we include some technical details on the treatment of the 
long-range interactions in our simulations. A preliminary 
account of part of this work has already appeared^. 



II. THEORETICAL APPROACH 

At each stage of the simulation the crystal is mod- 
eled as a slab of finite thickness. However, it is conve- 
nient to index the allowed B-sites as in an infinite three- 
dimensional crystal lattice 



I = i ai + j a 2 + k a 3 . 



(1) 



Two-dimensional (2-D) periodic boundary conditions 
(PBC) are employed along the a\ and a 2 direc- 
tions, which define the x-y Cartesian plane, using a 
L\a\ x £20*2 = Ai x A-2 2-D supercell. Growth pro- 
ceeds along the z-direction. The simulation is initialized 
as a slab of uniform thickness Hqcls, with a predefined 
B-atom configuration. A given simulation is terminated 
when cither the maximum slab thickness or the maximum 
number of Monte Carlo (MC) time steps is reached. We 
use the notation L x L x H max to label a particular simu- 
lation, where H max is the number of layers for maximum 
slab thickness (the initial substrate included). 

The SOS restriction that we impose does not allow 
the formation of voids. The crystal configuration, C, is 
specified at each stage of the simulation by the set of 



occupied sites I = (i,j,k) and their charges qi. The 
BV electrostatic model cannot be directly used in this 
slab geometry, due to ill-defined electrical boundary con- 
ditions in the z direction and the lack of exact charge 
neutrality during the growth simulation. Section A de- 
scribes how we handle these issues. Similarly, a direct 
application of the KMC algorithm is inefficient due to 
the long-range Coulomb interaction. Section B describes 
the KMC method and our modifications to make it ap- 
plicable to the model. 



A. The Electrostatic Model 

The BV model is derived by considering the total elec- 
trostatic energy for an A(BB'B )0 3 compound: 



(It,Vt') 



eRi T -Ri 



(2) 



where Ri T is the position of the ion on site t (=A, B, Oi, 
02, O3) in cell I, and e is the dielectric constant. For a 
given Bravais lattice, e sets the energy scale. We consider 
the perovskite structure with group II A-site atoms (e.g. 
Ba, Pb), so the charges on the A and O sites have fixed 
values of +2e and — 2e, respectively. Since the average B- 
site charge is +4e, it is convenient to express the charges 
on the B-sites, Qi,b, as 



11, B 



= 4e 



(3) 



Up to a constant, the configurationally averaged electro- 
static energy depends only on the B-site charges, since 
the configurational average of qi is zero: 



E B (C) = -J2 



(1,1' 



\i-v 



(4) 



where we have for simplicity restricted ourselves in 
Eq. Q to a cubic Bravais lattice with lattice parameter 
a, and Rib = la. In this model each cell I is therefore 
reduced to a single lattice site with charge qi , and the en- 
ergy of the compound is given by the inter-site Coulomb 
interaction. 

The long-range Coulomb interaction must be treated 
with care in a bulk simulation to ensure proper conver- 
gence. For 2-D and 3-D simulations with periodic bound- 
ary conditions, the method of Ewald is often used, in 
which periodic images of the charges and neutralizing 
background charges are introducedi^iSiiLiSiiS so that 
the bare Coulomb form — l'\ is replaced with a re- 
duced form v(l — I'). For our growth simulations, we are 
dealing with a slab geometry with PBC only in two di- 
mensions (x-y). Some modifications are required before 
the Ewald method can be applied. 

In the simulations, we will need to calculate the energy 
change from Eq. due to the evaporation of a charged 
ion qii at the surface of the crystal (see Eq. (|25|l below). 
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The distribution of point charges that qi' 
described by the charge density 



p(r)=££><$(r-*-fl) 



''sees" can be 



(5) 



l R 



where I runs through the position vectors of the atoms 
within the simulation cell, and R is a 2-D Bravais super- 
cell lattice vector: R = niAi + TI2A1. Directly summing 
the Coulomb potentials of the individual point charges, 
V(r — I — R) = qi/\r — I — R\, leads to an ill-defined 
and conditionally convergent result, as is well known. 
However, for three-dimensional periodic boundary con- 
ditions, a unique solution of Poisson's equation exists 
(for an electrically neutral system) , and it is conveniently 
calculated using Ewald's method. Subject to some addi- 
tional, physically motivated conditions, a unique solution 
can also be found for finite thickness slabs that are infi- 
nite in extent along two spatial directions. 

Solutions of Poisson's equation, W 2 V(r) = —A-Kp{r), 
in our simulations are subject to two-dimensional (2-D) 
PBC V(r+R) = V(r), as is the charge density p(r). The 
2-D PBC imply that V(r) and p{r) can be expanded as: 



P(r) = T,Pg{* 



„iG-r„ 



G 



V{r)=Y,V G {z)e i 

G 



(6) 



where G is a 2-D supcrcell reciprocal lattice vector and 
r p is the x-y component of r, r p = r—(r-z)z = iai+ja?. 
Substitution of Eqs. JBJ into Poisson's equation yields the 
ordinary differential equation 



d 2 V G (z) 
dz 2 



-G 2 V g (z) = -Ak Pg (z), 



whose solution can be expressed as 



00 

V G (z) = -4* J G(z - z')p G (z')dz\ 



(7) 



(8) 



where Q[z — z') is the Green's function corresponding to 
Eq. Q. 

If there are any ill-defined contributions to the 
Coulomb potential, they must arise from the G = so- 
lution in Eq. J7J) or (JSJ). This is because only the G = 
term of p(r) in Eqs. J5J contributes to the net slab charge. 
Even if the slab is electrically neutral, there may still 



be a net dipolc moment D, which would lead to differ- 
ent asymptotic values of Coulomb potential at z = ±00. 
Again, D also depends only on the G = term of p(r), 
where 



D 



zp(z)dz, 



(9) 



with 



p{z) = \j A p( r ) dxd y = pg=o(z), 



(10) 



where A is the area of the 2-D supercell. 

We therefore first consider the well-defined G ^ so- 
lutions of Eq. • Physically meaningful results require 
that the solutions satisfy lim^i^^ V G (z) — 0, which 
leads to the following unique definition of the G 7^ 
Green's function: 



Q(z-z') 



■d(z-z')e 



l\ p -G(z-z') 



+${z'-z)e G( - z - z,) 



2G 



(11) 



where G = \G\. For any reasonably localized charge dis- 
tribution Pg(z), Eqs. (JHJ) and result in well-behaved, 
exponentially decaying solutions of V G (z) as \z\ — > 00. 
For G = 0, Eq. Q becomes 



d 2 V (z) 
dz 2 



-4:TTp (z). 



(12) 



As adatoms are adsorbed or atoms evaporate in the 
course of the growth simulations, the net charge will 
fluctuate so that the total charge in the simulation su- 
percell will not be precisely zero at each stage of the 
simulation. Similarly a net dipolc D may form. How- 
ever, in a real growth process there are always compen- 
sating charges that will cancel any ill-defined long-range 
effects due to the lack of charge neutrality or the pres- 
ence of a dipolc moment. In our calculations, wc sim- 
ulate this by a construction that ensures that po(z) in 
Eq. I|12|) always represents a neutral charge distribution 
with D = 0. This leads to well-defined boundary condi- 
tions lim^i^oo V (z) = 0. 

As in the 3-D Ewald method, a diffuse localized charge 
density g(r) is added and subtracted to each point charge 
to facilitate the decomposition of the potential into ab- 
solutely convergent direct- and reciprocal lattice sums: 



p( r ) = EZ)«[*( r -'-- R )-5( r - , -- R )]+EZ)«f(r-i-ii) 

l R I R 

= pi(r)+p 2 (r). 

I 



(13) 



The diffuse charge density g(r) is chosen to be a normal- 



ized spherically symmetric Gaussian, as in the 3-D Ewald 
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method: 



V 7T 



3/2 



where the value of the Ewald convergence parameter a 
is arbitrary, but is usually chosen to optimize the con- 
vergence of both the direct- and reciprocal-lattice sums. 
The integrated charge of pi(r) is zero by construction, as 
is its dipolc moment D, so its contribution V\(r) to the 
Coulomb potential can be obtained by a rapidly conver- 
gent direct-lattice sum, given in the Appendix. 

The procedure for calculating the Coulomb potential 
V^{r) due to p2(f) requires special handling. V2(l'), the 
potential at the position of qi> in the simulation cell, is 
due to: i) the I ^ I' Gaussian charge densities and their 
periodic images qig(r — I — R), and ii) the periodic im- 
ages qi'g(r — I' — R). [As in the 3D Ewald method, a 
spurious interaction of the point charge qv with its own 
Gaussian density qi>g(i — V) is explicitly removed later.] 
Alternatively, the contribution (ii) above can be replaced 
by iia) the Gaussian densities —qig(r — V — R) located 
at the positions of the periodic images of V . In a bulk 
crystal simulation with 3-D PBC and a neutral simula- 
tion cell, these two formulations are equivalent, since the 
integrated total charge vanishes: 



where V denotes the projection of the position V onto 
the plane defined by the qi sublattice. The charge den- 

(14) sity p2 (r) has a rapidly convergent expansion in terms 
of 2-D planewaves given by Eq. 10. Moreover, the G = 

contribution of p 2 1 ' 1 ' (r) vanishes, so the Coulomb poten- 
tial V^(r) is readily found using Eqs. © and ifTTjl. These 
two approximations ensure overall average-charge neu- 
trality and vanishing dipolc moment D — 0, resulting 
in a well-defined Coulomb potential at each stage of the 
growth simulation. Complete formulas for the potential 
v(l'—l) are given in the Appendix. 



Ql' 



(15) 



In the 2-D slab geometry of our growth simulations, 
this will not be the case in general. Overall charge neu- 
trality is still satisfied in a statistical sense, however. Our 
procedure for calculating V^t*) consists of two approxi- 
mations. The first approximation is to use formulation 
(iia) above. This means that the contribution of each 
qig(r) sublattice to V2(V) is to be calculated as the po- 
tential due to the charge density: 



Ph 0) 



qi J2[9(r-l-R)-g(r-l'-R)}. (16) 

R 

Since the integrated charge of p 2 1 ' 1 \r) is zero, the use 

of this approximation effectively imposes overall charge 

neutrality at each stage of the growth simulation. 

The boundary conditions are still ill-defined however, 

(i i') 

since the sum of sublattice potentials due to the p\ ' (r) 
may still have a dipole moment D. We therefore in- 
troduce a second approximation: the Gaussian image 
densities —qig(r — I' — R) are made coplanar with the 
qig(r — I — R) sublattice. In other words, the Gaussian 
densities —qig(r) are placed at positions that are the pro- 
jections of the qi> image positions onto the plane defined 
by the qi sublattice. In place of Eq. (|16|) . the contribu- 
tion of each qig(r) sublattice is thus calculated as the 
potential due to the charge density: 



#°(r)=«£ 

R 



r-l-R)-g(r-l'-R) , (17) 



B. Kinetic Monte Carlo method for long-range 
interactions 

The kinetic Monte Carlo (KMC) method is one of 
several simulation techniques commonly employed to 
model the relaxation processes of systems away from 
equilibrium (e.g. growth processes). It has been ap- 
plied successfully to crystal growth and surface/interface 
phenomenalisms mostly in the context of kinetic Ising 
models. Due to the long-range interactions between ions 
in our electrostatic model, the usual implementation of 
KMC for Ising-likc models is inefficient, with the accep- 
tance rates of events becoming very low. We developed 
a modified sampling algorithm to make the simulation 
practical for this model. Here we briefly outline the basic 
theoretical background for the KMC method, and then 
describe our modifications and give the relevant imple- 
mentation details. 

In the KMC simulation, the dynamics of the system 
is described as stochastic processes such as adsorption, 
evaporation, and surface migration. We consider only the 
first two in our simulation. As mentioned, the adatoms 
represent the B-site ions in the single crystal perovskite 
alloy. They are characterized entirely by their charges 
and they interact with each other by the interaction de- 
scribed above. 

In the grand canonical ensemble, the Hamiltonian that 
will be used in the growth simulations can then be ex- 
pressed in term of Eq. Q as 



H(C) = E B (C) + ApN, 



(18) 



where N is the total number of adsorbed adatoms. The 
electrostatic energy term in the Hamiltonian is respon- 
sible for evaporation, while the second term, which de- 
pends on the chemical potential difference between the 
solid and the gas phases, controls the rate in which 
adatoms stick on the surface. The growth simulation 
is then characterized by competing adsorption and des- 
orption events. The SOS restriction imposed in the sim- 
ulation prevents formation of vacancies and allows us to 
write H as 



H(C) = E B (C) + AfiJ2k 



(19) 



ho 
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where hij is the number of layers in the present crystal 
configuration at the horizontal position ia\ + ja^. 

In KMC the time evolution of the system is simulated 
through a Markov chain of configurations. Let us de- 
fine P(C, t) time-dependent distribution of config- 
urations. The transition rate from C to C", a crystal 
configuration related to C by a single time step, is de- 
noted by w(C — > C). We then have the usual master 
equation^: 



dP(C, t) 
dt 



Y^w(C^C')P(C,t) 

c 

5>(C ->C)P(C,t), (20) 



C" 



where the first term on the right describes the loss be- 
cause of transitions away from C , while the second term 
describes the gain because of transitions into C. In the 
equilibrium limit (as t — > oo), the Boltzmann distribution 



Peg — Z 



exp 



-n(c) 

k R T 



(21) 



is reached, where k B is the Boltzmann constant. We 
require that detailed balance be satisfied: 



w{C^C) = P eq {C') 
w{C C) P eq (C) 



exp 



H(C')-H(C) 



k R T 



(22) 



The KMC technique can be viewed as a method of 
solving Eq. (|20fl stochastically. We adopt the following 
choice of transition rates w(C — ► C) 



W P 



exp (An/k B T) 
exp(-AE B (C)/k B T), 



(23) 
(24) 



where w a and w e are the rates for adsorption and evap- 
oration, respectively, of an adatom. It can be veri- 
fied straightforwardly that this choice indeed satisfies 
Eq. (|22|l . The rate w e for an adatom of charge q T i to 
evaporate from the surface depends on the change in to- 
tal potential energy in the crystal 



AE B (C) 



E B {C')-E B {C) 
qv_ 

ca 



(25) 



i 



For kinetic Ising models, the algorithm of BKLS allows 
an efficient stochastic realization of the kinetic process 
under the choice in Eq.'s l|23|l and l|24|l . In this algorithm, 
a site (i, j) is selected randomly in each step at the surface 
of the grown crystal. An event is then selected by Monte 
Carlo sampling^! from the list of three possible events, 
{adsorption, evaporation, nothing}. The interaction in 
Ising type models is limited to near-neighbors, and the 
energy difference AE B {C) is completely determined by 
the local environment at site The global maxi- 

mum of w e , i.e., the minimum possible energy change, 



A£ mm = mm[AE B (C)], can be obtained straightfor- 
wardly by considering all possible local configurations. 
This gives a corresponding global maximum of the evapo- 
ration rates: w™ ax = exp (— AE nun /k B T), which defines 
a normalization factor: 



W 



W a 



(26) 



The relative probabilities for the three events are there- 
fore 

w„ w P , . 

{P a = -^,P e = ^,P n = l-P a - Pe }. (27) 

With the electrostatic model, however, the energy 
change in Eq. I|24|) depends on the entire configuration 
C . It is therefore difficult to determine the global mini- 
mum, A_E mm . Indeed, even if AE mm could be identified, 
the energy change AE B (C), which can vary greatly with 
C and the simulation cell size, would be much greater 
than Ai? mm for most configurations. This would cause 
the evaporation and adsorption probabilities P a and P e 
to be small, with P n approaching unity. As a result the 
acceptance rate of events becomes small, and the algo- 
rithm becomes ineffective. 

To overcome this difficulty, we modify the standard 
algorithm so that all N — L\ x Li surface sites are con- 
sidered simultaneously, instead of sweeping through the 
surface sites. An event list is created which includes ev- 
ery possible event for every possible surface site. This in- 
creases the algorithm complexity, because of the need to 
store and update an array of surface potentials, calculate 
the event list, and sample an event from this list. The 
advantage is that an event is guaranteed to take place 
in each step of the algorithm and that the need for de- 
termining Ai? mm is eliminated. Evaporation/adsorption 
rates for all possible sites are normalized. The sum of 
the probabilities for an adsorption or evaporation to oc- 
cur at a surface site is unity. Specifically, the modified 
algorithm consists of the following steps: 

(i) Generate a list, £, of all possible events per time step. 
There are 2A^ possible events: an evaporation or an 
adsorption could happen on each of the N = L\ x L 2 
surface sites. 

(ii) Calculate the rates (w) of adsorption and evaporation 
for each site on the surface. Denote the total rates 

2N 

by W: W = J2m- 

i 

(iii) Normalize these 2N rates by W, giving probabil- 
ities, Pi, for adsorption and evaporation on sites 
1,2,---,2A. 

(iv) Generate a random number r G [0, 1) and choose the 

i 

first event E{ such that Pk > r - An event will 

fc=i 

always be chosen. 

(v) Generate the new configuration C based on the cho- 
sen event £ j. 
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(vi) Assign a "real time" increment Ai rea j = —\n(r')/W 
to this MC step, where r' is another random number 
on [0,1). 

The last step is a result of our considering the global 
event list and forcing an event to occur in every step. 
The issue of real "time" in a KMC simulation is a subtle 
one. Often the Monte Carlo time, £mc, is used as some 
measure of the real time. In the standard algorithm, the 
global normalization factor W (defined by u>™ ax ) controls 
the overall rate of events and sets a "time scale." In our 
approach, W is time-dependent, and an event is forced to 
happen in each step regardless of the total rate W for the 
configuration at hand. When W is low, an evaporation 
or adsorption is less likely to happen but one is selected 
anyway. Conversely, when W is high, an evaporation or 
adsorption is more likely to happen but still only one 
is selected. This introduces a bias which should vanish 
in the limit of large system size but which should be 
corrected for at finite L. Based on the rate equation, wc 
assume an exponential relation between time and W . A 
step in which W is high corresponds to a short time, and 
vice versa. Step (vi) is a way to account for this time scale 
stochastically, by rescaling AtMC with a MC sampling 
from an exponential distribution which is determined by 
the normalization factor W in each step. 

III. RESULTS 

We now present the results from our simulations for 
A(BB')03 and A(BB'b")0 3 crystals. Growth simula- 
tions are presented in Section A. Growth rates are stud- 
ied, and charge-charge correlation functions are calcu- 
lated to measure the degree of growth order. The effects 
of varying the crystallographic orientation of the slabs 
were explored, with the slabs labelled according to the 
slab perpendicular (z) direction. In A(BB B )03 sys- 
tems, a fraction of tetravalent B ions are mixed in. In 
our growth simulations, these tetravalent ions do not ap- 
pear to mix at low temperatures, choosing instead to 
phase-separate from the pure crystal. To further study 
this, Section B presents results of static total energy and 
free-energy calculations for fixed slab configurations. 

A. Crystal Growth 

The growth process is a function of temperature T, 
chemical potential difference A/i, and the Coulomb inter- 
action. These parameters are fixed throughout a given 
simulation. As discussed in Section Til Al and in the Ap- 
pendix, we tabulate v(V — I), and we will use reduced 
units in our simulations below. The energies ( A/z and 
Eb{C) ) are scaled by £ = 1/ea. There is only one free 
parameter between £ and the temperature fcgT, which 
sets the energy scale of the problem. Below we will give 
the temperature k B T in reduced units. For example, for 



a ~ 8 a.u. and e ~ 10 (typical values of BMN solid solu- 
tions) in Eq. J3J, 1350 C corresponds to k B T = 0.41 in 
the simulation. 

As an overview, Figs. ^ and [21 present a comparison of 
simulations of the simple III1/2V1/2 rocksalt alloy and a 
II 1 / 3 V 2 /3 hetcrovalent alloy such as BMN. (All substrates 
in our simulations have neutral surface layers.) We mea- 
sure the growth rate of the crystal based on the KMC 
dynamics. If Nq adatoms are gained in m MC steps 
(each defined as one attempt at the procedure outlined 
in Section lrTB|) . the growth rate is defined as 

T= v^m G A , 777- (28) 

w a 2^i=i Ai rca i(i) 

Note that as defined the growth rate T is renormalizcd 
by the absorption rate. The growth rate is plotted as a 
function of the chemical potential for a range of temper- 
atures. The rocksalt structure has layers of positive and 
negative charges alternating along the [111] direction. It 
typifies the crystal ordering of a wide variety of materials, 
including some of the perovskite alloys. Heterovalent bi- 
naries, described by II1/2VI1/2 (<?b = ±2) or IIIi/ 2 V 1 / 2 
(c[b = ±1), exhibit rocksalt B-site chemical order. By 
contrast, in the II1/3V2/3 heterovalent binary BMN the 
equilibrium state shows [lll]i : 2 ordering of two layers of 
metal group V(gs = +1) alternating with one layer of the 
group II(<7b = —2) atom. Both the rocksalt and BMN 
simulations were initialized with a 20-laycr thick slab, 
with perfect [lll]i : i and [lll]i : 2 ordering, respectively. 
The rocksalt simulation used a 2-D 12 x 12 supercell, 
while the BMN simulations were done mostly with 6x6 
supercells, although some simulations with 12 x 12 and 
15 x 15 were carried out to verify that the finite-size ef- 
fects were small. The rocksalt structure simulations ran 
for 1, 000L 2 MC steps, up to a maximum thickness of 100 
layers. For BMN, 10, 000L 2 MC steps were used, because 
for a given temperature and A^ growth was significantly 
slower. In Fig. [21 we show visualizations of the grown 
BMN crystals to illustrate the simulation environment 
and the 1:2 order at low temperatures with slow growth. 

The two sets of curves in Fig^andl^larc qualitatively 
similar. What is not evident from the figures, however, is 
the degree of order in each simulation. For a given tem- 
perature, as A/x increases, the adsorbtion rate in Eq. I|23|) 
increases, and adatoms are more likely to stick. For fixed 
Afi, as fcgT decreases, the adsorption rate will increase, 
but more importantly, the "selectiveness" of evaporation 
will increase. A lower k B T will, in effect, increase the en- 
ergy differences between competing configurations. The 
direct result, as growth is concerned, will be that adatoms 
will increasingly prefer to have more instead of less neigh- 
bors with correct charge ordering (layer-by-layer growth 
vs. rough growth), and adatoms with the same charge 
will seem more repulsive. For very high A/i, adatoms 
will stick anywhere, no matter what the location or ionic 
adversity is, and the growth rate will be high. Alterna- 
tively, if the temperature becomes too high, the crystal 
will melt, preferring the liquid phase, and result in neg- 
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FIG. 1: Rocksalt growth rate vs. chemical potential for a 
[001] slab. 
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FIG. 2: BMN growth rate vs. chemical potential for a [-1-11] 
slab. 



ative growth. 

To examine the degree of ordering, we computed the 
charge-charge correlation function. The Fourier trans- 
form of this correlation function, which we will denote 
by rj(k), gives the structure factor: 



V (k) 



a ^ llQl+l' cxp(-ifc • I') 
w 



(29) 



where a is the normalization factor, and k is the wave 
vector in the Brillouin zone of the unit cell. The magni- 
tude of rj(k) characterizes the B-site order, e.g., a large 
value of r\ at k = (|, |, i) indicates a strong [lll]i : i 
order while one at k = — 

a 

[111] i:2 order. 



i, i, i) indicates a strong 



(a) 



(b) 



FIG. 3: Visualizations of grown BMN crystals. Shown are 6 x 
6 supercells with: (a) growth direction along [001], fcsT = 0.1 
and A/i = —1.0; (b) growth direction along [111], fcsT = 0.1 
and A/i — —1.1. 



The growth rate and the charge-charge structure fac- 
tor in Eq. (gUJ are plotted in Figs. 14171 In each fig- 
ure, the displayed range of A/i was chosen to coincide 
with the range where the order parameter decreases from 
nearly unity (perfect order) to essentially zero (disorder) . 
As A/n increases the adsorption rate increases, but the 
growth is disordered and there is greater surface rough- 
ness. Indeed there is only a limited range where ordered 
growth occurs. The grown crystal structures are con- 
sistent with the observed ground state configuration of 
rocksalt (Fig. 0} and BMN (Fig. 's EE}. The most strik- 
ing difference between the growth behaviors of rocksalt 
and BMN is the enormous reduction of the growth rate of 
BMN compared to that of the rocksalt structure. More- 
over, for rocksalt the growth rate increases linearly as a 
function of A/i in the region where the order parameter r\ 
is rapidly decreasing. By contrast, the BMN growth rate 
is relatively constant in this region. As A/i increases be- 
yond this region, there is a sudden onset of much larger 
growth rates, but the resulting crystals are disordered. 
The growth rate of BMN increases as the temperature is 
increased f Figs. 15171 . 

We next attempted to model the growth of BMN-BZ 
(1-x) (Mg 1 / 3 Nb 2 /3) + x Zr solid solutions. In the elec- 
trostatic Hamiltonian in Eq. (0J, tetravalent Zr corre- 
sponds to a neutral charge qi = 0, so sites occupied 
by Zr have zero interaction energy. As in the simula- 
tions of pure BMN systems, the chemical composition 
determines the probabilities with which different charge 
species are adsorbed at the surface. In the initial sub- 
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FIG. 4: Rocksalt growth rate F of Eq. 1281 (top panel) and 1:1 
order parameter rj(k = —(5,5, §)) (bottom panel) vs. chem- 
ical potential. The temperature is ksT = 0.1 and the growth 
direction is [001]. A 12 x 12 supercell is used, with 1000 MC 
time steps. 
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ical potential. The temperature is ksT = 0.025 and the 
growth substrate direction is [111]. A 6 x 6 supercell is used, 
with 300,000 MC time steps. 
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FIG. 7: BMN growth rate and 1:2 order parameter 
vs. chemical potential. The temperature is fc_gT = 0.2. Other 
parameters are the same as in Fig. . 



strate, tetravalent ions with the corresponding concen- 
tration were incorporated, using random mixing (next 
section). With a l:2-ordered substrate, we studied con- 
centrations x ~ 10%, with temperatures of ksT ~ 0.1 
to 0.2, and varying the chemical potential A/i ~ —1.0 
to —0.5. Very little incorporation of the tetravalent ions 
occurred. We found similar results with an initially 1:1- 
ordered substrate (random-site model; sec below) , where 
a wider range of x was explored. Again the order of the 
substrate was not sufficient to induce the incorporation 
of tetravalent ions in the growth phase. Instead the sys- 
tem seemed to favor evaporating the adsorbed tetravalent 
ions more than the charged particles, to grow pure BMN. 



B. Energy Calculations 

To further study the inability to incorporate tetrava- 
lent ions at low temperatures, we examined the total en- 
ergy per particle En of fixed slab configurations of B-site 
order. A phase separated model, in which all the tetrava- 
lent adatoms were situated in the outermost surface lay- 
ers, was compared with various structural models that 
incorporated tetravalent ions. In each model, the calcu- 
lations were performed for two different configurational 
B-site orderings of the +2 and +5 ions (qi = — 2,+l, 
respectively). These configurations were the 1:1 and 1:2 
layering along [111] directions, i.e., [lll]i : i and [lll]i : 2 
order, repectively. 

The [lll]i:2 ordering corresponds to the x = order 
of BMN, with a layer of qi = —2 alternating with two 
layers of charge qi = +1 along the [111] direction. We 
chose the [lll]i : i ordering to correspond to the random- 
site model^, which is observed in the BMN-BZ equilib- 
rium simulations for x > 0.05m&2£ In the random-site 
model there are [111] layers of qi = +1 alternating with 
a mixed layer of charges qi = —2, +1, 0. The random-site 
model is meant to represent the presence of short-range 
B-site order from experimental observations. No long- 
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FIG. 8: Total energy per particle for B-site [lll]i : a ordering 
as a function of slab thickness 1/H and slab crystallographic 
orientation. Each set has three barely distinguishable curves, 
corresponding to three lattice sizes: 12 x 12, 15 x 15, and 
18 x 18. 



FIG. 9: Total energy per particle for [lll]i : 2 ordering, as a 
function of slab thickness H and (randomly mixed) tetrava- 
lent concentration x. E/N increases, as x increases from 0% 
to 2%, 5%, 10%, 15%, and 25%. A [111] slab with a 15 x 15 
supercell was used. 



range ordering has been observed. Nevertheless in our 
simple model here we will fix the ordered qi = +1 layers 
and choose the mixed layers to be a random mixture of 

We first examine finite-size effects in Fig.|S| which plots 
ejv as a function of slab thickness for various 2-D super- 
cells containing no tetravalent ions, for [111] 1: 2 order- 
ing. Results for [001] and [111] slabs are shown, both of 
which correspond to neutral surface layers. As H — > oo, 
£n ~ £n + const. /H as expected, where the constant 
e n represents the average bulk value and H is the slab 
thickness. 

Size effects with incorporated tetravalent ions are stud- 
ied next. Fig.Elplots en for [lll]i : 2 ordering as a function 
of slab thickness for various concentrations of randomly 
mixed tetravalent ions, using a [111] slab and 15 x 15 
supercell. These calculations are for a random distribu- 
tion of +0 (tetravalent) ions replacing -2 or +1 ions in 
an otherwise perfectly ordered [lll]i : 2 slab at each thick- 
ness H. Within statistical error bars, the asymptotic in- 
dependence is similar to that without tetravalent ions in 
Fig. |H1 Fig. EH plots en for [111] 1:1 ordering with and 
without randomly mixed 10% tetravalent ions. As seen 
in Fig- EI £n rapidly increases with increasing x. This 
is consistent with the inability to incorporate tetravalent 
ions in the growth simulations on [lll]i : 2 ordered slabs. 
By contrast en for [lll]x : i ordering is essentially inde- 
pendent of x within statistical error, as shown in Fig. 1101 

Fig.^Jplots En as a function of tetravalent concentra- 
tion x for random-mixing and phase-separation models, 
showing results for [lll]i : i and [lll]i : 2 ordered 12 x 12 
[001] slabs (H = 200). For the phase-separation model, 
the total number of ions includes the outermost layers 
of tetravalent ions. At x = the 1:2 ordered crystal 
has a lower energy than the 1:1 ordered crystal, which 
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FIG. 10: Total energy per particle vs. slab thickness. Results 
are shown for [lll]i : i ordering with (dashed line) and without 
(solid line) 10% randomly mixed tetravalent ions. A [001] slab 
with a 12 x 12 supercell was used. 



is consistent with our results from the growth simulation 
and with the observed ground state configuration of pure 
BMN. For random-mixing, e^v increases linearly with x 
for [111] i:2 ordering while it is essentially independent of 
x for [1 1 1] i:i ordering. In the phase-separation model, En 
increases linearly for both orderings. These results show 
that phase separation is favored for the [lll]i : 2 ordering, 
while random mixing is favored by [lll]i : i ordering. 

Fig. [n] illustrates why the growth simulations failed to 
incorporate tetravalent ions at low temperature. In the 
electrostatic model, the 1:2 ordered state is the ground 
state and is optimally ordered. The potential energy be- 
tween any charge and all other charges in the system is 
negative. For example, with a 18 x 18 slab this poten- 
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FIG. 11: Total energy per particle vs. tetravalent concentra- 
tion x for random-mixing and phase-separation models. Re- 
sults are shown for [lll]i : i and [lll]i : 2 ordered 12 x 12 [001] 
slabs (H = 200). 



tial energy is ~ —5.92 for a —2 charge and ~ —1.48 for 
a +1 charge. Thus, replacing a charge (either —2 or 
+1 ) by a neutral tetravalent ion in this state raises the 
total energy of the system, while a phase-separated con- 
figuration in which the tetravalent ion is placed away 
from the the ordered slab keeps the total energy un- 
changed. To examine this more closely, we calculated 
the free-energy (F = en — TS), where S is the mixing 
entropy due to the incorporated tetravalent ions. Fig.ll"2l 
plots the free-energy as a function of temperature for 
four concentrations of tetravalent ions. The free energy 
of the phase-separated 1:2 ordered slabs is constant in 
our model, because it is perfectly ordered and has van- 
ishing entropy. The free energy of the phase-separated 
1:1 ordered slabs decreases with increasing temperature, 
despite the perfectly ordered outmost layers of tetrava- 
lent ions, due to the mixing entropy of the random layers 
with —2, +1, and charges. In all cases in Fig. ^1 the 
phase-separated 1:2 ordered slabs have the lowest free en- 
ergy at low temperatures, where ordered crystal growth 
occurs in our simulations, but at temperatures between 
ksT ~ 1 — 2 the 1:2 ordered and the 1:1 ordered random 
mixing models start to be favored. 



IV. DISCUSSION 

There are striking differences between the growth be- 
havior of the III1/2V1/2 rocksalt ordered structure and 
the II1/3V2/3 BMN structure. The ordered rocksalt 
structure forms over a wide range of A/i (absorption 
rates) as shown in Fig. 4. By contrast, ordering of 
the 1:2 structure in BMN type crystals is more diffi- 
cult to achieve experimentallyiS*2i When these materials 
arc initially synthesized, they crystallize in a disordered 
structure. With extended annealing the 1:2 structure is 
approached t22> As discussed by Davies et alm&, the initial 
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FIG. 12: Free energy of BMN crystal for (a) 10% (b) 15% 
(c) 20% (d) 25% tetravalent concentrations. Symbols have 
the same meaning as in Fig. 1111 



synthesis and processing are controlled by irreversible ki- 
netic processes rather than by thermodynamic factors, 
and a more correct description of the formation of the 
1:2 ordered structures is in terms of the nucleation and 
growth of small ordered domains with increasing anneal- 
ing time and temperature. Eventually large (>100 nm) 
1:2 ordered domains are observed The need for long 
annealing times is consistent with our simulations. Figs. 
5-7 show that the range of A/i where ordered 1:2 growth 
occurs narrows as the temperature increases from UbT 
= 0.025 to 0.2. In this range, the growth rate is ap- 
proximately constant as a function of A/i. Moreover, 
when ordered crystal growth occurs, the BMN growth 
rate is much smaller than that of the rocksalt structure 
at the same temperature. Highly ordered growth was 
possible in the BMN simulations but required low tem- 
peratures and a delicate balance with the chemical po- 
tential. Neither of these requirements is likely to be met 
under experimental synthesis conditions. At tempera- 
tures corresponding to the actual sintering temperature 
of BMN (kgT ~ 0.5), large growth rates can be achieved, 
as shown in Fig. 2, but the growth is highly disordered. 
The long annealing times allow the slow formation of 
the 1:2 ordered regions. In our KMC simulations, diffu- 
sion processes are excluded so there can be no annealing. 
We also note that the growth rate was sensitive to the 
slab orientation. For example, we found that growth rate 
along [111] direction was almost an order of magnitude 
larger than that along [001]. 

Our results arc also qualitatively consistent with the 
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long experimental history of failed attempts to coarsen 
the 1:1 ordered nanoscale domains in PMN type crys- 
tals. Prior to the experiments of Akbas and Davies^, the 
1:1 ordered regions were apparently limited to nanoscale 
size and represented only a small volume fraction of the 
crystal. The space-charge model, which was invoked 
to explain this behavior, hypothesized that the 1:1 or- 
dered regions arose from a rocksalt ordering of the -2 
and +1 B-site charges, implying charge-imbalanced 1:1 
domains. The apparently limited size of these domains 
could be explained by the rapidly increasing energy of 
larger domains due to Coulomb repulsion. With care- 
ful annealling at much higher temperatures than had 
previously been tried, however, some fully 1:1 ordered 
crystals were synthesized 2 ^. Our calculations show that 
long-range ionic interactions favor the growth of disor- 
dered crystals, and ordering occurs only after annealing. 
Moreover, ionic interactions appear to favor the 1:2 or- 
dering. However, entropic contributions to the free en- 
ergy and short-range covalent interactions tend to favor 
1:1 ordering. Covalent bonding is negligible for Ba ions 
but very important for Pb ions. Thus there is a delicate 
competition between 1:2 and 1:1 ordering for doping with 
small concentrations of the tetravalent ions in (l-x)BMN- 
xBZ and (l-x)PMN-xPT. In (l-x)BMN-xBZ, there is a 
crossover from 1:2 to 1:1 ordering as x increases to about 
5%. While in (l-x)PMN-xPT, the stronger short-range 
covalent bonding of Pb favors 1:1 ordering at all concen- 
trations. 

For pure systems, our minimal paradigm for growth 
simulations captures the differences in growth rate and 
ordering between rocksalt-type and BMN-type crystal 
growth. This indicates that the simple ionic model is 
a reasonable starting point for describing the growth of 
perovskite solid solutions. More direct and quantitative 
comparisons with experiment will require additional in- 
gredients such as short-range interactions and the inclu- 
sion of diffusive processes. 

For systems with tetravalent ions, our results show that 
the ground state is a phase-separated state of tetravalent 
ions and 1:2 ordered BMN over a wide range of tetrava- 
lent compositions. On the other hand, equilibrium sim- 
ulations of the ionic mode l 10 ' 22 suggest that for x > 0.05 
the 1:1 ordering is preferred, with no phase separation. 
Several factors distinguish these calculations, which likely 
have to do with the apparent contradiction in their ob- 
servations. The first is the difference in the nature of 
the simulations. In our growth simulation, tetravalent 
ions are allowed to evaporate from the crystal, which fa- 
cilitates phase separation. The equilibrium calculations 
were done in the canonical ensemble with the tetravalent 
ions mixed in, where it is more difficult to detect phase 
separation without large simulation cell sizes. Our simu- 
lations were at lower temperatures where ordered growth 
could be induced by tuning the chemical potential Afi 
(absorption rate). At these temperatures the system is 
essentially in the ground state, as Fig.[3|shows. Incorpo- 
ration of tetravalent ions could be induced at larger A/x, 



which is expected as adsorption dominates evaporation, 
but in this case random growth occurs. Secondly, since 
our [lll]i:i structure is an artificial model of random mix- 
ing of —2, +1, and neutral charges in one layer and per- 
fectly ordered +1 in another, its energy must be higher 
than the actual 1:1 structure achieved in the equilibrium 
simulations. This means that the actual cross-over of 
the random-mixing [lll]i : i structure will occur at lower 
temperatures. Indeed, the fcgT ~ 0.25 equilibrium calcu- 
lations show [lll]i : i ordering for concentrations x greater 
than about 0.05. Thus the absence of phase separation 
in the equilibrium calculations might be due to a lower 
free-energy than our estimate in Fig. 1121 from the arti- 
ficial random-site structure. Our results combined with 
the equilibrium calculations therefore suggest the follow- 
ing picture of the equilibrium state of the ionic model. In 
the ground state phase-separation takes place for x > 0. 
Beyond some x-dependent critical temperature tetrava- 
lent ions are incorporated, most likely in a structure that 
favors 1:1 order. 

To determine if the new phase (phase-separation) at 
low temperatures that we have found is realistic for these 
alloys, the ionic model must be improved. One possibil- 
ity is first-principles based H c g, which have shown great 
promise in describing ferroelectrics and simple solid- 
solutions^. Like the Ising model these H c g project out 
what are considered to be the most important ionic de- 
grees of freedom. In addition to the long-range Coulomb 
interaction, short-range interactions are also included. 
The H e g parameters are fitted to the results of a set of 
first-principles density-functional calculations, so there is 
effectively no experimental input (except sometimes the 
average crystal volume). The simplified form of i? e ff for 
ferroelectrics and ferroelectric alloys has permitted sim- 
ulations of equilibrium properties on thousands of atoms 
as a function of temperature and applied external elec- 
tric field. A main difficulty in applying these in a growth 
simulation is computational cost, which has typically re- 
quired fixed distributions of B-site ions even in equilib- 
rium simulations of solid-solutions. In our kinetic Monte 
Carlo model, another possibly important factor that is 
not included is surface diffusion. Coupled with the solid- 
on-solid restriction, the simulation is severely limited in 
its ability to "heal" disorder, and these approximations 
may have contributed to low ordered growth rates and 
raised the critical temperature for phase separation. Re- 
moval of these restrictions would improve the model and 
increase its applicability. 



V. SUMMARY 

The growth of the technologically important BMN 
type perovskite alloys was studied by kinetic Monte Carlo 
using an ionic model. An enhanced KMC algorithm was 
formulated to treat long-range Coulomb interactions effi- 
ciently. We found that this minimal paradigm was capa- 
ble of describing ordering features of the growth of pure 
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BMN and PMN type single crystals. The largest growth 
rates were observed along the [III] direction, but best 
ordered growth rates are substantially less than those 
of rocksalt. Highly ordered growth was possible, but re- 
quired very low temperatures and a delicate balance with 
the chemical potential. For mixed systems such as BMN- 
BZ, we found that the T = ground state of the model 
was one in which tetravalent ions phase separate from 
a 1:2 ordered pure system. As a result, little incorpo- 
ration of tetravalent ions occurs in the growth process 
at low temperatures. At higher temperatures, tetrava- 
lent ions can be incorporated, but the resulting crystals 
show no chemical ordering. The tendency of the purely 
ionic model to favor phase separate was further stud- 
ied using free energy calculations determined from T = 
total energy calculations and including a mixing en- 
tropy. This indicated that, if diffusive mechanisms were 
included, chemical orderings consistent with those found 
in equilibrium studies could develop at the higher tem- 
peratures characteristic of realistic alloy synthesis. 
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APPENDIX A: COULOMB POTENTIAL 

The 2-D Ewald potential, is given as the sum of three 
terms 



V(l' -l)= Vl(V V 2 {i - l) + V S (1'), 



(Al) 



where v\ and v 2 are due to pi(r) and p 2 ix), respectively 
in Eq. I|13|) . and v s is the correction for the interaction 
of the point charge qi> with its own Gaussian density 
q v g(r - V) in p 2 (r). 



To calculate Vi(l' — I) we place, for consistency, the 
(R 7^ 0) qi> images at their vertical projections onto the 
plane of the qi sublattice. V\{V — I) is then given by 



a/ ,\ erfc( v / a|Z'— I — R 
= \V-l-R\ 
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crfc 




I'-l-R 
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(A2) 



The mathematical form of this contribution is identical 
to its 3-D counterpart, except that the sum is over 2-D 
rather than 3-D direct-lattice vectors R. 

The 2-D planewave expansion of p 2 1 (r) in Eq. I|17|) 
is given by 



Pi ( r ) = 1i 



G#0 



-G 2 /4a 



-iGla 



1/2 



-iG.l' 



iG- 



(A3) 



where we have used the fact that l' z = l z . Substituting 
into Eq. JHJ and using Eq. (fTT|) . yields: 



v 2 (i'-i)= £ -m[f(G)-f{-G)] 



e iG-(l' p -l p )_ 1 



where 



(A4) 



(A5) 



Finally, the correction for the interaction of the point 
charge qi> with its own Gaussian density is given by: 



erf(y / a 


V - 


-V) 




V - 


-V 





(A6) 



As verified by direct calculation, the sum of these three 
terms in independent of the parameter a. For efficiency, 
v{V — I) is stored look-up table. 
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